rm(list=ls())
library(dplyr) # v.1.0.7
library(psy) # v.1.1
library(ggplot2) # v.3.3.5
library(stargazer) # v.5.2.2
library(factoextra) # v.1.0.7
library(polycor) # v.0.8.0
library(mediation) # v.4.5.0
library(psych) # v.2.1.9
library(GPArotation) # v.2014.11.1
library(lavaan) # v.0.6.9
library(nFactors) # v.2.4.1
library(sjPlot) # v.2.8.10
library(sjmisc) # v.2.8.9
library(ggExtra)  # v.0.9
library(interplot) # v.0.2.3
library(gridExtra) # v.2.3
library(jtools) # v.2.1.4
library(dotwhisker) # v.0.7.4
library(ggpubr) # v.0.4.0

############################################################################
# Yellow Peril or Model Minority? Measuring Janus-Faced Prejudice toward Asians in the United States
# D.G. Kim & Enze Han
# Political Science Research and Methods
# 2024
############################################################################

############################################################################
# Setup & data loading
## Set working directory to your file where the following two files are saved
lucid <- read.csv("Lucid_clean.csv")
dynata <- read.csv("Dynata_clean.csv")

############################################################################
# Saving entire output of main replication file to a log

sink("replication_log.txt", split = TRUE)  # "split = TRUE" will display output in the console as well
source("Kim & Han PSRM_final.R") # sink() is placed at the end to close the loop
sink()

############################################################################
# Replication codes below are for the figures and tables in the main paper first (Figures 2-6 & Table 2)
# followed by those in the online appendix
############################################################################

############################################################################
# Starting with the figures & tables in the main paper
############################################################################

############################################################################
# Figure 2: Confirmatory factor analysis and inter-item correlations

## Left panel in Figure 2: confirmatory factor analysis 
## (using Dynata sample; results with the Lucid sample are in the appendix and replicated below)

aarmms.model <- 'f1 =~ aar1 + aar2 + aar3 + aar4
f2 =~ mms1 + mms2 + mms3 + mms4'

fit.dynata <- cfa(aarmms.model, data=dynata, std.lv = TRUE,ordered = FALSE)
summary(fit.dynata, fit.measures=T, standardized=T)
fitmeasures(fit.dynata, c('cfi', 'rmsea', 'rmsea.pvalue', 'bic', 'chisq', 'df', 'srmr'))  
parameterEstimates(fit.dynata, standardized=TRUE, rsquare=TRUE) # standardized factor loadings are under "est" column with rows of "op = =~"

## Right panel in Figure 2: Pearson correlation coefficients among the eight (AAR & MMS) scale items
## (using Dynata sample; results with the Lucid sample are in the appendix and replicated below)

library(GGally)
p <- ggpairs(dynata, columns = c("aar1", "aar2", "aar3", "aar4", "mms1", "mms2", "mms3", "mms4"), 
        title = "", switch="both",
        columnLabels = c('AAR_1', 'AAR_2', 'AAR_3', 'AAR_4', 'MMS_1', 'MMS_2', 'MMS_3', 'MMS_4'),
        diag = list(continuous = wrap("densityDiag")),
        lower = list(continuous = wrap("cor", color="black", size=6)),                                                                                                                                                                            
        upper = list(continuous = wrap("blank",alpha = 0.3,size = 1, color="black"))) + theme_bw()

ggsave("figure 2.png", plot = p)

############################################################################
# Figure 3: Distributions of and correlations between AAR and MMS

## Left panel in Figure 3: scatter & density plots (using Dynata sample)
s2 <- ggscatterhist(dynata, x = "mms", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"), margin.params = list(fill = "lightgray"),
                    margin.plot=c("density"),
                    margin.plot.size = 1.5, margin.ggtheme = theme_minimal(),
                    add = "loess", conf.int = TRUE,
                    cor.coef = TRUE, cor.method = "pearson",color="#000000",
                    xlab = "MMS", ylab = "AAR", cor.coef.size = 6,
                    title="Dynata sample") + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

ggsave("figure 3_left.png", plot = s2)

## Right panel in Figure 3: distributions of responses to the two scales split by the midpoint with raw counts in parentheses

s4 <- ggscatterhist(dynata, x = "mms", y = "aar", size=.01,point=FALSE, cor.coeff.args = list(label.sep = "\n"), margin.params = list(fill = "lightgray"),
                    margin.plot=c("density"),
                    margin.plot.size = 1.5, margin.ggtheme = theme_minimal(),
                    add = "none", conf.int = FALSE,
                    cor.coef = FALSE, cor.method = "pearson",color="#000000",
                    xlab = "MMS", ylab = "AAR", cor.coef.size = 6,
                    title="Dynata sample") + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

ggsave("figure 3_right.png", plot = s4)

prop.table(table(dynata$aar>3 & dynata$mms>3)) # 25.6% (n=259) are in the upper right corner
table(dynata$aar>3 & dynata$mms>3) # 25.6% (n=259) are in the upper right corner

prop.table(table(dynata$aar>3 & dynata$mms<=3)) # 5.4% (n=55) are in the upper left corner
table(dynata$aar>3 & dynata$mms<=3) # 5.4% (n=55) are in the upper left corner

prop.table(table(dynata$aar<=3 & dynata$mms<=3)) # 27.9% (n=282) are in the lower left corner
table(dynata$aar<=3 & dynata$mms<=3) # 27.9% (n=282) are in the lower left corner

prop.table(table(dynata$aar<=3 & dynata$mms>3)) # 41.0% (n=414) are in the lower right corner
table(dynata$aar<=3 & dynata$mms>3) # 41.0% (n=414) are in the lower right corner

############################################################################
# Figure 4: Demographic and dispositional correlates of AAR (Dynata sample)

dynata$stereo <- (dynata$violent_asian + (8-dynata$lazy_asian) + (8-dynata$intel_asian) + (8-dynata$trust_asian))*(1/4)
dynata$stereo_white <- (dynata$violent_white + (8-dynata$lazy_white) + (8-dynata$intel_white) + (8-dynata$trust_white))*(1/4)
dynata$stereo_black <- (dynata$violent_black + (8-dynata$lazy_black) + (8-dynata$intel_black) + (8-dynata$trust_black))*(1/4)
dynata$stereo_hispanic <- (dynata$violent_hispanic + (8-dynata$lazy_hispanic) + (8-dynata$intel_hispanic) + (8-dynata$trust_hispanic))*(1/4)

library(ggpubr)

c1 <- ggscatter(dynata, x = "feeling_4", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,
                cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                xlab = "Favorability", ylab = "AAR", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c2 <- ggscatter(dynata, x = "stereo", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,
                cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                xlab = "Negative stereotype", ylab = "AAR", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))


c3 <- ggscatter(dynata, x = "ethno", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,
                cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                xlab = "Ethnocentrism", ylab = "AAR", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))


c4 <- ggscatter(dynata, x = "sdo", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,fill = "lightgrey",
                cor.coef = TRUE, cor.method = "pearson",
                xlab = "SDO", ylab = "AAR", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))


c5 <- ggscatter(dynata, x = "whiteidentity", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,
                cor.coef = TRUE, cor.method = "pearson", fill = "lightgrey",
                xlab = "White identity", ylab = "AAR", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))


c6 <- ggscatter(dynata, x = "finance1", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,
                cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                xlab = "Financial stress", ylab = "AAR", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))


c7 <- ggscatter(dynata, x = "ideo", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,
                cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                xlab = "Ideology", ylab = "AAR", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))


c8 <- ggscatter(dynata, x = "pid", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,
                cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                xlab = "Party ID", ylab = "AAR", cor.coef.size = 6)+ ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))


c9 <- ggscatter(dynata, x = "sr", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,
                cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                xlab = "Symbolic racism", ylab = "AAR", cor.coef.size = 6)+ ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))


c10 <- ggscatter(dynata, x = "edu", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "Education", ylab = "AAR", cor.coef.size = 6)+ ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))


c11 <- ggscatter(dynata, x = "age", y = "aar", size=.1, point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "Age", ylab = "AAR", cor.coef.size = 6)+ ylim(1,5) + xlim(20,80) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c12 <- ggscatter(dynata, x = "inc", y = "aar", size=.01, point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "Income", ylab = "AAR", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

ggsave("figure_4.png", plot = arrangeGrob(c1, c2, c3, c4, c5, c9, c7, c8, c6, c10, c11, c12, ncol=4), width = 10, height = 8, dpi = 300)


############################################################################
# Figure 5: Demographic and dispositional correlates of MMS

c13 <- ggscatter(dynata, x = "feeling_4", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "Favorability", ylab = "MMS", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c14 <- ggscatter(dynata, x = "stereo", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "Negative stereotype", ylab = "MMS", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c15 <- ggscatter(dynata, x = "ethno", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "Ethnocentrism", ylab = "MMS", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c16 <- ggscatter(dynata, x = "sdo", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "SDO", ylab = "MMS", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c17 <- ggscatter(dynata, x = "whiteidentity", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "White identity", ylab = "MMS", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c18 <- ggscatter(dynata, x = "finance1", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "Financial stress", ylab = "MMS", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c19 <- ggscatter(dynata, x = "ideo", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "Ideology", ylab = "MMS", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))


c20 <- ggscatter(dynata, x = "pid", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "Party ID", ylab = "MMS", cor.coef.size = 6)+ ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c21 <- ggscatter(dynata, x = "sr", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "Symbolic racism", ylab = "MMS", cor.coef.size = 6)+ ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c22 <- ggscatter(dynata, x = "edu", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "Education", ylab = "MMS", cor.coef.size = 6)+ ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c23 <- ggscatter(dynata, x = "age", y = "mms", size=.1, point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "Age", ylab = "MMS", cor.coef.size = 6)+ ylim(1,5) + xlim(20,80) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c24 <- ggscatter(dynata, x = "inc", y = "mms", size=.01, point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",
                 xlab = "Income", ylab = "MMS", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c25 <- grid.arrange(c13, c14, c15, c16, c17, c21, c19, c20, c18, c22, c23, c24, ncol=4)

ggsave("figure_5.png", plot = arrangeGrob(c13, c14, c15, c16, c17, c21, c19, c20, c18, c22, c23, c24, ncol=4), width = 10, height = 8, dpi = 300)


############################################################################
# Figure 6: Cognitive responses to AAR and MMS scale items
## Tables A.3 and A.4 in the appendix are replicated here, since Figure 6 is derived from estimates in those two tables.
### First, run the codes in the descriptive statistics section and find replications below

### Descriptive statistics
#### Coder 1
table(dynata$aar_prej1)
prop.table(table(dynata$aar_prej1))

table(dynata$mms_prej1)
prop.table(table(dynata$mms_prej1))

table(dynata$aar_neg1)
prop.table(table(dynata$aar_neg1))

table(dynata$mms1_neg1)
prop.table(table(dynata$mms_neg1))

table(dynata$aar_pos1)
prop.table(table(dynata$aar_pos1))

table(dynata$mms_pos1)
prop.table(table(dynata$mms_pos1))

#### Coder 2
table(dynata$aar_prej2)
prop.table(table(dynata$aar_prej2))

table(dynata$mms_prej2)
prop.table(table(dynata$mms_prej2))

table(dynata$aar_neg2)
prop.table(table(dynata$aar_neg2))

table(dynata$mms1_neg2)
prop.table(table(dynata$mms_neg2))

table(dynata$aar_pos2)
prop.table(table(dynata$aar_pos2))

table(dynata$mms_pos2)
prop.table(table(dynata$mms_pos2))

#### Aggregate
table(dynata$aar_prej_b)
prop.table(table(dynata$aar_prej_b)) # ~37% mention the statement as prejudiced at least once for AAR

table(dynata$mms_prej_b)
prop.table(table(dynata$mms_prej_b)) # ~49% mention the statement as prejudiced at least once for MMS

table(dynata$aar_neg_b)
prop.table(table(dynata$aar_neg_b)) # ~15% mention negative traits at least once for AAR

table(dynata$mms_neg_b)
prop.table(table(dynata$mms_neg_b)) # ~5% mention negative traits at least once for MMS

table(dynata$aar_pos_b)
prop.table(table(dynata$aar_pos_b)) # ~39% mention positive traits at least once for AAR

table(dynata$mms_pos_b)
prop.table(table(dynata$mms_pos_b)) # ~39% mention positive traits at least once for MMS

table(dynata$aar_indi_b)
prop.table(table(dynata$aar_indi_b)) # ~33% affirm individualism for AAR

table(dynata$mms_indi_b)
prop.table(table(dynata$mms_indi_b)) # ~35% affirm individualism for MMS

#### Regressions (Probit)
dynata$pid2 <- factor(dynata$pid2)
dynata$race <- factor(dynata$race)
dynata$male <- factor(dynata$male)

lucid$pid2 <- factor(lucid$pid2)
lucid$race <- factor(lucid$race)
lucid$male <- factor(lucid$male)

summary(o1 <- glm(aar_prej_b ~ aar + pid2 + ideo + edu + male + race + age + inc, family = binomial(link = "probit"), data = dynata))
summary(o2 <- glm(aar_neg_b ~ aar + pid2 + ideo + edu + male + race + age + inc, family = binomial(link = "probit"), data = dynata))
summary(o3 <- glm(aar_pos_b ~ aar + pid2 + ideo + edu + male + race + age + inc, family = binomial(link = "probit"), data = dynata))
summary(o4 <- glm(aar_indi_b ~ aar + pid2 + ideo + edu + male + race + age + inc, family = binomial(link = "probit"), data = dynata))
summary(oo4 <- glm(aar_indi_b ~ sr + pid2 + ideo + edu + male + race + age + inc, family = binomial(link = "probit"), data = dynata))

summary(o5 <- glm(mms_prej_b ~ mms + pid2 + ideo + edu + male + race + age + inc, family = binomial(link = "probit"), data = dynata))
summary(o6 <- glm(mms_neg_b ~ mms + pid2 + ideo + edu + male + race + age + inc, family = binomial(link = "probit"), data = dynata))
summary(o7 <- glm(mms_pos_b ~ mms + pid2 + ideo + edu + male + race + age + inc, family = binomial(link = "probit"), data = dynata))
summary(o8 <- glm(mms_indi_b ~ mms + pid2 + ideo + edu + male + race + age + inc, family = binomial(link = "probit"), data = dynata))
summary(oo8 <- glm(mms_indi_b ~ sr + pid2 + ideo + edu + male + race + age + inc, family = binomial(link = "probit"), data = dynata))

summary(o9 <- glm(aar_prej_b ~ sdo + pid2 + ideo + edu + male + race + age + inc, family = binomial(link = "probit"), data = dynata))
summary(o10 <- glm(mms_prej_b ~ sdo + pid2 + ideo + edu + male + race + age + inc, family = binomial(link = "probit"), data = dynata))
summary(o11 <- glm(aar_prej_b ~ whiteidentity + pid2 + ideo + edu + male + age + inc, family = binomial(link = "probit"), data = dynata))
summary(o12 <- glm(mms_prej_b ~ whiteidentity + pid2 + ideo + edu + male + age + inc, family = binomial(link = "probit"), data = dynata))

#### Table A.3: Predictors of key themes in open-ended responses to AAR scale items from Dynata sample
stargazer(type="text", o1,o2,o3,o4,oo4,o9,o11,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

#### Table A.4: Predictors of key themes in open-ended responses to MMS scale items from Dynata sample
stargazer(type="text", o5,o6,o7,o8,oo8, o10,o12,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

#### Figure 6 (main paper): Cognitive responses to AAR and MMS scale items

po1 <- plot_model(o1, type = "pred", terms = c("aar"), show.data=FALSE, colors="bw", ci.lvl=.95) + ylim(0,1) +
  labs(title="Prejudiced statement", y ="") +
  xlab("AAR") +
  theme(title=element_text(size=16),
        axis.title.x = element_text(size = 18), 
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(color="black", size=14, angle=0),
        axis.text.x = element_text(color="black", size=12, angle=0))

po2 <- plot_model(o2, type = "pred", terms = c("aar"), show.data=FALSE, colors="bw", ci.lvl=.95) + ylim(0,1) +
  labs(title="Negative trait mentioned", y ="") +
  xlab("AAR") +
  theme(title=element_text(size=16),
        axis.title.x = element_text(size = 18), 
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(color="black", size=14, angle=0),
        axis.text.x = element_text(color="black", size=12, angle=0))

po3 <- plot_model(o3, type = "pred", terms = c("aar"), show.data=FALSE, colors="bw", ci.lvl=.95) + ylim(0,1) +
  labs(title="Positive trait mentioned", y ="") +
  xlab("AAR") +
  theme(title=element_text(size=16),
        axis.title.x = element_text(size = 18), 
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(color="black", size=14, angle=0),
        axis.text.x = element_text(color="black", size=12, angle=0))

po4 <- plot_model(o4, type = "pred", terms = c("aar"), show.data=FALSE, colors="bw", ci.lvl=.95) + ylim(0,1) +
  labs(title="Individualism affirmed", y ="") +
  xlab("AAR") +
  theme(title=element_text(size=16),
        axis.title.x = element_text(size = 18), 
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(color="black", size=14, angle=0),
        axis.text.x = element_text(color="black", size=12, angle=0))

poo4 <- plot_model(oo4, type = "pred", terms = c("sr"), show.data=FALSE, colors="bw", ci.lvl=.95) + ylim(0,1) +
  labs(title="Individualism affirmed", y ="") +
  xlab("Symbolic racism") +
  theme(title=element_text(size=16),
        axis.title.x = element_text(size = 18), 
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(color="black", size=14, angle=0),
        axis.text.x = element_text(color="black", size=12, angle=0))

po5 <- plot_model(o5, type = "pred", terms = c("mms"), show.data=FALSE, colors="bw", ci.lvl=.95) + ylim(0,1) +
  labs(title="Prejudiced statement", y ="") +
  xlab("MMS") +
  theme(title=element_text(size=16),
        axis.title.x = element_text(size = 18), 
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(color="black", size=14, angle=0),
        axis.text.x = element_text(color="black", size=12, angle=0))

po6 <- plot_model(o6, type = "pred", terms = c("mms"), show.data=FALSE, colors="bw", ci.lvl=.95) + ylim(0,1) +
  labs(title="Negative trait mentioned", y ="") +
  xlab("MMS") +
  theme(title=element_text(size=16),
        axis.title.x = element_text(size = 18), 
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(color="black", size=14, angle=0),
        axis.text.x = element_text(color="black", size=12, angle=0))

po7 <- plot_model(o7, type = "pred", terms = c("mms"), show.data=FALSE, colors="bw", ci.lvl=.95) + ylim(0,1) +
  labs(title="Positive trait mentioned", y ="") +
  xlab("MMS") +
  theme(title=element_text(size=16),
        axis.title.x = element_text(size = 18), 
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(color="black", size=14, angle=0),
        axis.text.x = element_text(color="black", size=12, angle=0))

po8 <- plot_model(o8, type = "pred", terms = c("mms"), show.data=FALSE, colors="bw", ci.lvl=.95) + ylim(0,1) +
  labs(title="Individualism affirmed", y ="") +
  xlab("MMS") +
  theme(title=element_text(size=16),
        axis.title.x = element_text(size = 18), 
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(color="black", size=14, angle=0),
        axis.text.x = element_text(color="black", size=12, angle=0))

poo8 <- plot_model(oo8, type = "pred", terms = c("sr"), show.data=FALSE, colors="bw", ci.lvl=.95) + ylim(0,1) +
  labs(title="Individualism affirmed", y ="") +
  xlab("Symbolic racism") +
  theme(title=element_text(size=16),
        axis.title.x = element_text(size = 18), 
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(color="black", size=14, angle=0),
        axis.text.x = element_text(color="black", size=12, angle=0))

po9 <- plot_model(o9, type = "pred", terms = c("sdo"), show.data=FALSE, colors="bw", ci.lvl=.95) + ylim(0,1) +
  labs(title="Prejudiced statement", y ="") +
  xlab("SDO") +
  theme(title=element_text(size=16),
        axis.title.x = element_text(size = 18), 
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(color="black", size=14, angle=0),
        axis.text.x = element_text(color="black", size=12, angle=0))

po10 <- plot_model(o10, type = "pred", terms = c("sdo"), show.data=FALSE, colors="bw", ci.lvl=.95) + ylim(0,1) +
  labs(title="Prejudiced statement", y ="") +
  xlab("SDO") +
  theme(title=element_text(size=16),
        axis.title.x = element_text(size = 18), 
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(color="black", size=14, angle=0),
        axis.text.x = element_text(color="black", size=12, angle=0))

po11 <- plot_model(o11, type = "pred", terms = c("whiteidentity"), show.data=FALSE, colors="bw", ci.lvl=.95) + ylim(0,1) +
  labs(title="Prejudiced statement", y ="") +
  xlab("White identity importance") +
  theme(title=element_text(size=16),
        axis.title.x = element_text(size = 18), 
        axis.title.y = element_text(size = 18),
        axis.text.y = element_text(color="black", size=14, angle=0),
        axis.text.x = element_text(color="black", size=12, angle=0))

po12 <- plot_model(o12, type = "pred", terms = c("whiteidentity"), show.data=FALSE, colors="bw", ci.lvl=.95) + ylim(0,1) +
  labs(title="Prejudiced statement", y ="") +
  xlab("White identity importance") +
  theme(
    title=element_text(size=16),
    axis.title.x = element_text(size = 18), 
    axis.title.y = element_text(size = 18),
    axis.text.y = element_text(color="black", size=14, angle=0),
    axis.text.x = element_text(color="black", size=12, angle=0))

blankPlot <- ggplot()+geom_blank(aes(1,1))+
  theme(
    plot.background = element_blank(), 
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(), 
    panel.border = element_blank(),
    panel.background = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(), 
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    axis.line = element_blank()
  )

#### Figure 6 Left panel
ggsave("figure_6_left.png", plot = arrangeGrob(po2, po3, po1, blankPlot,po9, po11, po4, poo4, ncol=2), width = 10, height = 8, dpi = 300)

#### Figure 6 Right panel
ggsave("figure_6_right.png", plot = arrangeGrob(po6, po7, po5, blankPlot, po10, po12, po8, poo8, ncol=2), width = 10, height = 8, dpi = 300)


############################################################################
# Table 2: AAR, MMS, and support for pro-Asian public policies

# Set up ###################################################################
## Standardizing key predictors first
dynata$aar_z <- scale(dynata$aar, center = TRUE, scale = TRUE) 
lucid$aar_z <- scale(lucid$aar, center=TRUE, scale=TRUE)
dynata$mms_z <- scale(dynata$mms, center = TRUE, scale = TRUE)
lucid$mms_z <- scale(lucid$mms, center=TRUE, scale=TRUE)
dynata$ideo_z <- scale(dynata$ideo, center = TRUE, scale = TRUE)
lucid$age_z <- scale(lucid$age, center=TRUE, scale=TRUE)
dynata$age_z <- scale(dynata$age, center=TRUE, scale=TRUE)
lucid$inc_z <- scale(lucid$inc, center=TRUE, scale=TRUE)
dynata$inc_z <- scale(dynata$inc, center=TRUE, scale=TRUE)
lucid$edu_z <- scale(lucid$edu, center=TRUE, scale=TRUE)
dynata$edu_z <- scale(dynata$edu, center=TRUE, scale=TRUE)

## Coding key outcome variables
### Immigrants from Asia
lucid$immi <- 1-lucid$asianpolicy_1
dynata$immi <- 1-dynata$asianpolicy_1

### Asian Americans elected in public offices
lucid$elec <- 1-lucid$asianpolicy_2
dynata$elec <- 1-dynata$asianpolicy_2

### Asian American students admitted to top U.S. universities
lucid$univ <- 1-lucid$asianpolicy_4
dynata$univ <- 1-dynata$asianpolicy_4

### Federal spending on programs on anti-Asian hate crimes
lucid$hate <- 1-lucid$asianpolicy_3

### Chinese students studying in the U.S.
dynata$chin <- 1-dynata$asianpolicy_3
############################################################################

## Now let's delve into Table 2. This table is composed of multiple, separate regression results:
## "Each alone," 
## "Both together," 
## "AAR x MMS" (coding details below)

### Each alone: AAR only (Full results shown in Table A.5.1 in the appendix)
summary(rr1 <- lm(immi ~ aar_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr2 <- lm(immi ~ aar_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr3 <- lm(elec ~ aar_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr4 <- lm(elec ~ aar_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr5 <- lm(univ ~ aar_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr6 <- lm(univ ~ aar_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr7 <- lm(hate ~ aar_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr8 <- lm(chin ~ aar_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

### Each alone: AAR only (Table A.5.2: controlling for symbolic racism and fiscal conservatism))
summary(rr1 <- lm(immi ~ aar_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr2 <- lm(immi ~ aar_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr3 <- lm(elec ~ aar_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr4 <- lm(elec ~ aar_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr5 <- lm(univ ~ aar_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr6 <- lm(univ ~ aar_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr7 <- lm(hate ~ aar_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr8 <- lm(chin ~ aar_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

### Each alone: MMS only (Full results shown in Table A.6.1 in the appendix)
summary(rr1 <- lm(immi ~ mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr2 <- lm(immi ~ mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr3 <- lm(elec ~ mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr4 <- lm(elec ~ mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr5 <- lm(univ ~ mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr6 <- lm(univ ~ mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr7 <- lm(hate ~ mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr8 <- lm(chin ~ mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

### Each alone: MMS only (Table A.6.2: controlling for symbolic racism and fiscal conservatism))
summary(rr1 <- lm(immi ~ mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr2 <- lm(immi ~ mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr3 <- lm(elec ~ mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr4 <- lm(elec ~ mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr5 <- lm(univ ~ mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr6 <- lm(univ ~ mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr7 <- lm(hate ~ mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr8 <- lm(chin ~ mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))


### Both together (Full results shown in Table A.7.1 in the appendix)
summary(rr1 <- lm(immi ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr2 <- lm(immi ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr3 <- lm(elec ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr4 <- lm(elec ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr5 <- lm(univ ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr6 <- lm(univ ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr7 <- lm(hate ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr8 <- lm(chin ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

### Both together (Table A.7.2: controlling for symbolic racism and fiscal conservatism))
summary(rr1 <- lm(immi ~ aar_z + mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr2 <- lm(immi ~ aar_z + mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr3 <- lm(elec ~ aar_z + mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr4 <- lm(elec ~ aar_z + mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr5 <- lm(univ ~ aar_z + mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr6 <- lm(univ ~ aar_z + mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr7 <- lm(hate ~ aar_z + mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr8 <- lm(chin ~ aar_z + mms_z + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

### AAR x MMS  (Full results shown in Table A.8.1 in the appendix)

#### Setup first: coding combinations of [high vs. low AAR] x [high vs. low MMS] ("strongly agree & agree" = high)
dynata$aar_high[dynata$aar>3] <- 1
dynata$aar_high[dynata$aar<=3] <- 0
table(dynata$aar_high)

dynata$mms_high[dynata$mms>3] <- 1
dynata$mms_high[dynata$mms<=3] <- 0
table(dynata$mms_high)

dynata$aarmms[dynata$aar_high==0 & dynata$mms_high==0] <- 0
dynata$aarmms[dynata$aar_high==0 & dynata$mms_high==1] <- 1
dynata$aarmms[dynata$aar_high==1 & dynata$mms_high==0] <- 2
dynata$aarmms[dynata$aar_high==1 & dynata$mms_high==1] <- 3
dynata$aarmms <- factor(dynata$aarmms)
table(dynata$aarmms)

lucid$aar_high[lucid$aar>3] <- 1
lucid$aar_high[lucid$aar<=3] <- 0
table(lucid$aar_high)

lucid$mms_high[lucid$mms>3] <- 1
lucid$mms_high[lucid$mms<=3] <- 0
table(lucid$mms_high)

lucid$aarmms[lucid$aar_high==0 & lucid$mms_high==0] <- 0
lucid$aarmms[lucid$aar_high==0 & lucid$mms_high==1] <- 1
lucid$aarmms[lucid$aar_high==1 & lucid$mms_high==0] <- 2
lucid$aarmms[lucid$aar_high==1 & lucid$mms_high==1] <- 3
lucid$aarmms <- factor(lucid$aarmms)
table(lucid$aarmms)


library(gmodels)
CrossTable(dynata$aar_high, dynata$mms_high)
CrossTable(lucid$aar_high, lucid$mms_high)

#### Full results shown in Table A.8.1 in the appendix

summary(rr1 <- lm(immi ~ aarmms + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr2 <- lm(immi ~ aarmms + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr3 <- lm(elec ~ aarmms + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr4 <- lm(elec ~ aarmms + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr5 <- lm(univ ~ aarmms + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr6 <- lm(univ ~ aarmms + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr7 <- lm(hate ~ aarmms + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr8 <- lm(chin ~ aarmms + sr_z  + smallgovernment_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001)) 

#### Table A.8.2 in the appendix

summary(rr1 <- lm(immi ~ aarmms + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr2 <- lm(immi ~ aarmms + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr3 <- lm(elec ~ aarmms + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr4 <- lm(elec ~ aarmms + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr5 <- lm(univ ~ aarmms + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr6 <- lm(univ ~ aarmms + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr7 <- lm(hate ~ aarmms + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr8 <- lm(chin ~ aarmms + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001)) 


############################################################################
#######     ! Now the figures & tables in the online appendix !    #########
############################################################################

############################################################################
# Figure A.1.1: Confirmatory factor analysis and correlation matrix (Lucid)

## Left panel in Figure A.1.1: confirmatory factor analysis 
## (using Lucid sample)

lucid$aar1 <- 6-lucid$aar_1 # reverse-coding AAR and MMS items
lucid$aar2 <- 6-lucid$aar_2
lucid$aar3 <- 6-lucid$aar_3
lucid$aar4 <- 6-lucid$aar_4
lucid$mms1 <- 6-lucid$mms_1
lucid$mms2 <- 6-lucid$mms_2
lucid$mms3 <- 6-lucid$mms_3
lucid$mms4 <- 6-lucid$mms_4

aarmms.model <- 'f1 =~ aar1 + aar2 + aar3 + aar4
f2 =~ mms1 + mms2 + mms3 + mms4'

fit.lucid <- cfa(aarmms.model, data=lucid, std.lv = TRUE,ordered = FALSE)
summary(fit.lucid, fit.measures=T, standardized=T)
fitmeasures(fit.lucid, c('cfi', 'rmsea', 'rmsea.pvalue', 'bic', 'chisq', 'df', 'srmr', 'tli'))  
parameterEstimates(fit.lucid, standardized=TRUE, rsquare=TRUE) # standardized factor loadings are under "est" column with rows of "op = =~"

## Right panel in Figure A.1.1: Pearson correlation coefficients among the eight (AAR & MMS) scale items
## (using Lucid sample)

library(GGally)
p <- ggpairs(lucid, columns = c("aar1", "aar2", "aar3", "aar4", "mms1", "mms2", "mms3", "mms4"), 
        title = "", switch="both",
        columnLabels = c('AAR_1', 'AAR_2', 'AAR_3', 'AAR_4', 'MMS_1', 'MMS_2', 'MMS_3', 'MMS_4'),
        diag = list(continuous = wrap("densityDiag")),
        lower = list(continuous = wrap("cor", color="black", size=6)),                                                                                                                                                                            
        upper = list(continuous = wrap("blank",alpha = 0.3,size = 1, color="black"))) + theme_bw()

ggsave("figure a.1.1.png", plot = p)

############################################################################
# Figure A.1.2: Principal component analysis and exploratory factor analysis

AARMMS_lucid <-with(lucid, data.frame(aar_1, aar_2, aar_3, aar_4, mms_1, mms_2, mms_3, mms_4))
AARMMS_dynata <-with(dynata, data.frame(aar_1, aar_2, aar_3, aar_4, mms_1, mms_2, mms_3, mms_4))

scree1 <- scree(AARMMS_lucid)
parallel1 <- fa.parallel(AARMMS_lucid, ylabel="Eigen Values of Principal Components and Factors")
parallel2 <- fa.parallel(AARMMS_dynata, ylabel="Eigen Values of Principal Components and Factors")

## Saving image outputs
png("figure a.1.2_left.png", width = 800, height = 600)
fa.parallel(AARMMS_lucid, ylabel = "Eigen Values of Principal Components and Factors")
dev.off()

png("figure a.1.2_right.png", width = 800, height = 600)
fa.parallel(AARMMS_dynata, ylabel = "Eigen Values of Principal Components and Factors")
dev.off()

############################################################################
# Figure A.2: Distributions of and correlations between AAR and MMS (Lucid sample)

## Left panel in Figure A.2: scatter & density plots (using Lucid sample)
s2 <- ggscatterhist(lucid, x = "mms", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"), margin.params = list(fill = "lightgray"),
                    margin.plot=c("density"),
                    margin.plot.size = 1.5, margin.ggtheme = theme_minimal(),
                    add = "loess", conf.int = TRUE,
                    cor.coef = TRUE, cor.method = "pearson",color="#000000",
                    xlab = "MMS", ylab = "AAR", cor.coef.size = 6,
                    title="Dynata sample") + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

ggsave("figure a.2_left.png", plot = s2)

## Right panel in Figure A.2: distributions of responses to the two scales split by the midpoint with raw counts in parentheses

s4 <- ggscatterhist(lucid, x = "mms", y = "aar", size=.01,point=FALSE, cor.coeff.args = list(label.sep = "\n"), margin.params = list(fill = "lightgray"),
                    margin.plot=c("density"),
                    margin.plot.size = 1.5, margin.ggtheme = theme_minimal(),
                    add = "none", conf.int = FALSE,
                    cor.coef = FALSE, cor.method = "pearson",color="#000000",
                    xlab = "MMS", ylab = "AAR", cor.coef.size = 6,
                    title="Dynata sample") + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

ggsave("figure a.2_right.png", plot = s4)

prop.table(table(lucid$aar>3 & lucid$mms>3)) # 32.8% (n=606) are in the upper right corner
table(lucid$aar>3 & lucid$mms>3) # 32.8% (n=606) are in the upper right corner

prop.table(table(lucid$aar>3 & lucid$mms<=3)) # 7.6% (n=141) are in the upper left corner
table(lucid$aar>3 & lucid$mms<=3) # 7.6% (n=141) are in the upper left corner

prop.table(table(lucid$aar<=3 & lucid$mms<=3)) # 24.4% (n=451) are in the lower left corner
table(lucid$aar<=3 & lucid$mms<=3) # 24.4% (n=451) are in the lower left corner

prop.table(table(lucid$aar<=3 & lucid$mms>3)) # 35.1% (n=649) are in the lower right corner
table(lucid$aar<=3 & lucid$mms>3) # 35.1% (n=649) are in the lower right corner


############################################################################
# Figure A.3: Demographic and dispositional correlates of AAR (Lucid sample)

library(ggpubr)
c1 <- ggscatter(lucid, x = "feeling_4", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,
                cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey", color="#000000",
                xlab = "Favorability", ylab = "AAR", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c4 <- ggscatter(lucid, x = "sdo", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,fill = "lightgrey",color="#000000",
                cor.coef = TRUE, cor.method = "pearson",
                xlab = "SDO", ylab = "AAR", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c5 <- ggscatter(lucid, x = "whiteidentity", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,
                cor.coef = TRUE, cor.method = "pearson", fill = "lightgrey",color="#000000",
                xlab = "White identity", ylab = "AAR", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c6 <- ggscatter(lucid, x = "finance1", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,
                cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                xlab = "Financial stress", ylab = "AAR", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c7 <- ggscatter(lucid, x = "ideo", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,
                cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                xlab = "Ideology", ylab = "AAR", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))


c8 <- ggscatter(lucid, x = "pid", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,
                cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                xlab = "Party ID", ylab = "AAR", cor.coef.size = 6)+ ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c9 <- ggscatter(lucid, x = "sr", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                add = "loess", conf.int = TRUE,
                cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                xlab = "Symbolic racism", ylab = "AAR", cor.coef.size = 6)+ ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c10 <- ggscatter(lucid, x = "edu", y = "aar", size=.01,point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                 xlab = "Education", ylab = "AAR", cor.coef.size = 6)+ ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c11 <- ggscatter(lucid, x = "age", y = "aar", size=.1, point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                 xlab = "Age", ylab = "AAR", cor.coef.size = 6)+ ylim(1,5) + xlim(1,6) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c12 <- ggscatter(lucid, x = "inc", y = "aar", size=.01, point=TRUE, cor.coeff.args = list(label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                 xlab = "Income", ylab = "AAR", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

grid.arrange(c1, c4, c5, c9, c7, c8, c6, c10, c11, c12, ncol=4)
ggsave("figure_a.3.png", plot = arrangeGrob(c1, c4, c5, c9, c7, c8, c6, c10, c11, c12, ncol=4), width = 10, height = 8, dpi = 300)

############################################################################
# Figure A.4: Demographic and dispositional correlates of MMS (Lucid sample)

c13 <- ggscatter(lucid, x = "feeling_4", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                 xlab = "Favorability", ylab = "MMS", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c16 <- ggscatter(lucid, x = "sdo", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                 xlab = "SDO", ylab = "MMS", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c17 <- ggscatter(lucid, x = "whiteidentity", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                 xlab = "White identity", ylab = "MMS", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c18 <- ggscatter(lucid, x = "finance1", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                 xlab = "Financial stress", ylab = "MMS", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c19 <- ggscatter(lucid, x = "ideo", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                 xlab = "Ideology", ylab = "MMS", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))


c20 <- ggscatter(lucid, x = "pid", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                 xlab = "Party ID", ylab = "MMS", cor.coef.size = 6)+ ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c21 <- ggscatter(lucid, x = "sr", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                 xlab = "Symbolic racism", ylab = "MMS", cor.coef.size = 6)+ ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c22 <- ggscatter(lucid, x = "edu", y = "mms", size=.01,point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                 xlab = "Education", ylab = "MMS", cor.coef.size = 6)+ ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c23 <- ggscatter(lucid, x = "age", y = "mms", size=.1, point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                 xlab = "Age", ylab = "MMS", cor.coef.size = 6)+ ylim(1,5) + xlim(1,6) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

c24 <- ggscatter(lucid, x = "inc", y = "mms", size=.01, point=TRUE, cor.coeff.args = list(label.y = 1.5,label.sep = "\n"),
                 add = "loess", conf.int = TRUE,
                 cor.coef = TRUE, cor.method = "pearson",fill = "lightgrey",color="#000000",
                 xlab = "Income", ylab = "MMS", cor.coef.size = 6) + ylim(1,5) +
  theme(axis.title.x = element_text(size = 16), axis.title.y = element_text(size = 16),
        axis.text.x = element_text(color="black", size=11, angle=0), 
        axis.text.y = element_text(color="black", size=11, angle=0))

grid.arrange(c13, c16, c17, c21, c19, c20, c18, c22, c23, c24, ncol=4)
ggsave("figure_a.3.png", plot = arrangeGrob(c13, c16, c17, c21, c19, c20, c18, c22, c23, c24, ncol=4), width = 10, height = 8, dpi = 300)


############################################################################
# Figure A.5: Demographic and dispositional correlates of MMS (Lucid sample)

dynata$sr1 <- 6-dynata$sr_1
dynata$sr2 <- dynata$sr_2
dynata$sr3 <- dynata$sr_3
dynata$sr4 <- 6-dynata$sr_4

lucid$sr1 <- 6-lucid$sr_1
lucid$sr2 <- lucid$sr_2
lucid$sr3 <- lucid$sr_3
lucid$sr4 <- 6-lucid$sr_4

library(GGally)
## Upper panel in Figure A.5 (Dynata sample)
p <- ggpairs(dynata, columns = c("mms1", "mms2", "mms3", "mms4", "sr1", "sr2", "sr3", "sr4"), 
        title = "Dynata sample", 
        columnLabels = c('MMS_1', 'MMS_2', 'MMS_3', 'MMS_4', 'SR_1', 'SR_2', 'SR_3', 'SR_4'),
        diag = list(continuous = wrap("densityDiag")),
        upper = list(continuous = wrap("cor", color="black", size=4)),                                                                                                                                                                            
        lower = list(continuous = wrap("smooth",alpha = 0.01,size = 1, color="black", lwd=0.1, method="loess", se=T))) + theme_bw()

ggsave("figure a.5_upper.png", plot = p)

## Lower panel in Figure A.5 (Lucid sample)
p <- ggpairs(lucid, columns = c("mms1", "mms2", "mms3", "mms4", "sr1", "sr2", "sr3", "sr4"), 
        title = "Lucid sample", 
        columnLabels = c('MMS_1', 'MMS_2', 'MMS_3', 'MMS_4', 'SR_1', 'SR_2', 'SR_3', 'SR_4'),
        diag = list(continuous = wrap("densityDiag")),
        upper = list(continuous = wrap("cor", color="black", size=4)),                                                                                                                                                                            
        lower = list(continuous = wrap("smooth",alpha = 0.01,size = 1, color="black", lwd=0.1, method="loess", se=T))) + theme_bw()

ggsave("figure a.5_lower.png", plot = p)

############################################################################
# Table A.9.1: Effects of feeling thermometer and stereotype measures

## Setup: coding relative scores for feeling thermometer
lucid$diff01 <- lucid$feeling_4 - lucid$feeling_1 # Asian - White
lucid$diff02 <- lucid$feeling_4 - lucid$feeling_2 # Asian - Black
lucid$diff03 <- lucid$feeling_4 - lucid$feeling_3 # Asian - Hispanic

lucid$feeling <- ifelse(lucid$race=="1", lucid$diff01,
                        ifelse(lucid$race=="2", lucid$diff02,
                               ifelse(lucid$race=="3", lucid$diff03, NA)))

lucid$feeling_z <- scale(lucid$feeling, center=TRUE, scale=TRUE)

dynata$diff01 <- dynata$feeling_4 - dynata$feeling_1 # Asian - White
dynata$diff02 <- dynata$feeling_4 - dynata$feeling_2 # Asian - Black
dynata$diff03 <- dynata$feeling_4 - dynata$feeling_3 # Asian - Hispanic

dynata$feeling <- ifelse(dynata$race=="1", dynata$diff01,
                         ifelse(dynata$race=="2", dynata$diff02,
                                ifelse(dynata$race=="3", dynata$diff03, NA)))
dynata$feeling_z <- scale(dynata$feeling, center=TRUE, scale=TRUE)

summary(rr1 <- lm(immi ~ feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr2 <- lm(immi ~ feeling_z + violent_z + untrust_z + intel_z + hard_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr3 <- lm(elec ~ feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr4 <- lm(elec ~ feeling_z + violent_z + untrust_z + intel_z + hard_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr5 <- lm(univ ~ feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr6 <- lm(univ ~ feeling_z + violent_z + untrust_z + intel_z + hard_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr7 <- lm(hate ~ feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr8 <- lm(chin ~ feeling_z + violent_z + untrust_z + intel_z + hard_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))

## This replicates estimates in Table A.9.1
stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001)) 

############################################################################
# Table A.9.2: Effects of AAR, MMS, and feeling thermometer I

summary(rr1 <- lm(immi ~ aar_z + mms_z + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr2 <- lm(immi ~ aar_z + mms_z + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr3 <- lm(elec ~ aar_z + mms_z + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr4 <- lm(elec ~ aar_z + mms_z + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr5 <- lm(univ ~ aar_z + mms_z + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr6 <- lm(univ ~ aar_z + mms_z + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr7 <- lm(hate ~ aar_z + mms_z + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr8 <- lm(chin ~ aar_z + mms_z + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))

## This replicates estimates in Table A.9.2
stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

############################################################################
# Table A.9.3: Effects of AAR, MMS, and feeling thermometer II

summary(rr1 <- lm(immi ~ aarmms + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr2 <- lm(immi ~ aarmms + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr3 <- lm(elec ~ aarmms + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr4 <- lm(elec ~ aarmms + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr5 <- lm(univ ~ aarmms + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr6 <- lm(univ ~ aarmms + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))
summary(rr7 <- lm(hate ~ aarmms + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=lucid))
summary(rr8 <- lm(chin ~ aarmms + feeling_z + pid2 + ideo_z + edu + male  + age + inc + race, data=dynata))

## This replicates estimates in Table A.9.3
stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

############################################################################
# Table A.10: AAR, MMS, and support for African American racial policies

dynata$blackpolicy_1 <- 1-dynata$blackpolicy_1 
dynata$blackpolicy_2 <- 1-dynata$blackpolicy_2
dynata$blackpolicy_3 <- 1-dynata$blackpolicy_3 
lucid$blackpolicy_1 <- 1-lucid$blackpolicy_1 
lucid$blackpolicy_2 <- 1-lucid$blackpolicy_2
lucid$blackpolicy_3 <- 1-lucid$blackpolicy_3 

summary(br01 <- lm(blackpolicy_1 ~ sr_z + aar_z +mms_z + smallgovernment_z + pid2 + ideo_z + edu + male  + race + age + inc, data=lucid))
summary(br02 <- lm(blackpolicy_1 ~ sr_z + aar_z +mms_z + smallgovernment_z + pid2 + ideo_z + edu + male  + race + age + inc, data=dynata))
summary(br03 <- lm(blackpolicy_2 ~ sr_z + aar_z +mms_z + smallgovernment_z + pid2 + ideo_z + edu + male  + race + age + inc, data=lucid))
summary(br04 <- lm(blackpolicy_2 ~ sr_z + aar_z +mms_z + smallgovernment_z + pid2 + ideo_z + edu + male  + race + age + inc, data=dynata))
summary(br05 <- lm(blackpolicy_3 ~ sr_z + aar_z +mms_z + smallgovernment_z + pid2 + ideo_z + edu + male  + race + age + inc, data=lucid))
summary(br06 <- lm(blackpolicy_3 ~ sr_z + aar_z +mms_z + smallgovernment_z + pid2 + ideo_z + edu + male  + race + age + inc, data=dynata))

stargazer(type="text", br01, br02, br03, br04, br05, br06,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

############################################################################
# Table A.11: Self-monitoring scale and social desirability bias

lucid$monitor <- (lucid$monitor1 + lucid$monitor2 + lucid$monitor3)*(1/3)
lucid$monitor_r <- (lucid$monitor-1)*(1/4)

lucid$aar_r <- (lucid$aar-1)*(1/4)
lucid$mms_r <- (lucid$mms-1)*(1/4)

summary(sm1 <- lm(mms_r ~ monitor + pid2 + ideo_z + edu + male  + race + age + inc, data=lucid))
summary(sm2 <- lm(aar_r ~ monitor + pid2 + ideo_z + edu + male  + race + age + inc, data=lucid))
summary(sm3 <- lm(monitor_r ~ factor(aarmms) + pid2 + ideo_z + edu + male  + race + age + inc, data=lucid))

stargazer(type="text", sm1, sm2, sm3,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

############################################################################
# Table A.14.1: White vs. non-white subgroup analysis I

## Dynata sample
t.test(dynata$aar[dynata$race==1], dynata$aar[dynata$race==3]) # White vs. Hispanic
t.test(dynata$aar[dynata$race==1], dynata$aar[dynata$race==2]) # White vs. Black

## Lucid sample
t.test(lucid$aar[lucid$race==1], lucid$aar[lucid$race==3]) # White vs. Hispanic
t.test(lucid$aar[lucid$race==1], lucid$aar[lucid$race==2]) # White vs. Black

############################################################################
# Table A.14.2: White vs. non-white subgroup analysis II (Table 2 with Blacks only)

### Each alone: AAR only
summary(rr1 <- lm(immi ~ aar_z + pid2 + ideo_z + edu + male  + age + inc, data=lucid, subset=race==2))
summary(rr2 <- lm(immi ~ aar_z + pid2 + ideo_z + edu + male  + age + inc, data=dynata, subset=race==2))
summary(rr3 <- lm(elec ~ aar_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr4 <- lm(elec ~ aar_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))
summary(rr5 <- lm(univ ~ aar_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr6 <- lm(univ ~ aar_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))
summary(rr7 <- lm(hate ~ aar_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr8 <- lm(chin ~ aar_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

### Each alone: MMS only
summary(rr1 <- lm(immi ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr2 <- lm(immi ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))
summary(rr3 <- lm(elec ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr4 <- lm(elec ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))
summary(rr5 <- lm(univ ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr6 <- lm(univ ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))
summary(rr7 <- lm(hate ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr8 <- lm(chin ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

### Both together
summary(rr1 <- lm(immi ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr2 <- lm(immi ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))
summary(rr3 <- lm(elec ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr4 <- lm(elec ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))
summary(rr5 <- lm(univ ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr6 <- lm(univ ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))
summary(rr7 <- lm(hate ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr8 <- lm(chin ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

### AAR x MMS 

summary(rr1 <- lm(immi ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr2 <- lm(immi ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))
summary(rr3 <- lm(elec ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr4 <- lm(elec ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))
summary(rr5 <- lm(univ ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr6 <- lm(univ ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))
summary(rr7 <- lm(hate ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==2))
summary(rr8 <- lm(chin ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==2))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))
############################################################################
# Table A.14.3: White vs. non-white subgroup analysis III (Table 2 with Hispanics only)

### Each alone: AAR only
summary(rr1 <- lm(immi ~ aar_z + pid2 + ideo_z + edu + male  + age + inc, data=lucid, subset=race==3))
summary(rr2 <- lm(immi ~ aar_z + pid2 + ideo_z + edu + male  + age + inc, data=dynata, subset=race==3))
summary(rr3 <- lm(elec ~ aar_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr4 <- lm(elec ~ aar_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))
summary(rr5 <- lm(univ ~ aar_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr6 <- lm(univ ~ aar_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))
summary(rr7 <- lm(hate ~ aar_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr8 <- lm(chin ~ aar_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

### Each alone: MMS only
summary(rr1 <- lm(immi ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr2 <- lm(immi ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))
summary(rr3 <- lm(elec ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr4 <- lm(elec ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))
summary(rr5 <- lm(univ ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr6 <- lm(univ ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))
summary(rr7 <- lm(hate ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr8 <- lm(chin ~ mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

### Both together
summary(rr1 <- lm(immi ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr2 <- lm(immi ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))
summary(rr3 <- lm(elec ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr4 <- lm(elec ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))
summary(rr5 <- lm(univ ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr6 <- lm(univ ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))
summary(rr7 <- lm(hate ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr8 <- lm(chin ~ aar_z + mms_z + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001))

### AAR x MMS 

summary(rr1 <- lm(immi ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr2 <- lm(immi ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))
summary(rr3 <- lm(elec ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr4 <- lm(elec ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))
summary(rr5 <- lm(univ ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr6 <- lm(univ ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))
summary(rr7 <- lm(hate ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=lucid, subset=race==3))
summary(rr8 <- lm(chin ~ aarmms + pid2 + ideo_z + edu + male  + age + inc , data=dynata, subset=race==3))

stargazer(type="text", rr1,rr2,rr3,rr4, rr5,rr6,rr7,rr8,
          title="",
          dep.var.labels =c(),
          keep.stat = c("n", "adj.rsq"), 
          digits=2, star.cutoffs = c(0.05, 0.01, 0.001)) 


sink()

############################################################################
#######                         Thank you                          #########
############################################################################